#!/usr/bin/python -itt

######### standard library imports ###########################################
import struct
import sys
import time 
from pprint import pprint
import logging
import logging.config
import os

######### external library imports ###########################################
import numpy as np
import matplotlib.pyplot as plt
from ROOT import *

######### project library imports ############################################
import new_helpers as helpers


scriptpath = os.path.dirname(os.path.realpath(__file__))
logging.config.fileConfig(scriptpath+'/logging.conf')
logger = logging.getLogger('reflector')


def findMirrorIntersection( run, dictname='photonData', 
            pos_key='ground_pos_turned', dir_key='ground_dir_turned',
            mirror_key='mirrors',
            new_pos_key='pos_on_mirror', new_mirror_id_key='mirror_id', new_time_key='time_to_mirror'):
    """ find intersection points of photons with mirrors.
    
    AS input we need:
    
    * the photon position   --> *pos_key* = 'ground_pos_turned'
    * the photon direction  --> *dir_key* = 'ground_dir_turned'
    * a list of mirrors     --> *mirror_key* = 'mirrors'
    
    Also new keys will be inserted into the dict given by *dictname*. The newly inserted info will be:
    
    * the position of the photon on the mirror --> *new_pos_key* = 'pos_on_mirror'
    * the ID of this mirror (an int) --> *new_mirror_id_key* = 'mirror_id'
    * the time, the photon needed from ground to the mirror --> *new_time_key* = 'time_to_mirror'
    
    The photons, that did not hit any mirror are being deleted.
    """
    
    start = time.time()
    logger.debug("findMirrorIntersection starting up ----------")
    data = run[dictname]
    mirrors = run['settings'][mirror_key]
    
    # we are going to add some information to the run.
    # this will be:
    #   * where a photon hit a mirror
    #   * which mirror it was
    #   * flighttime between mirror and ground (in direction of photon)
    #
    # so lets prepare the run for it
    
    
    num_photons = data[pos_key].shape[0]
    
    logger.debug("   Adding new keys to run['%s']: %s" % (dictname,str([new_pos_key, new_mirror_id_key, new_time_key]) ) )
    
    data[new_pos_key] = np.zeros( (num_photons,3), dtype=np.float32)
    data[new_mirror_id_key] = np.empty( num_photons, dtype=np.float32)
    data[new_mirror_id_key][:] = np.NaN
    data[new_time_key] = np.zeros( num_photons, dtype=np.float64)
    
    # the line defined by the photon is used to find the intersection point 
    # with the plane of each facette mirror. Then I check,
    # if the intersection point lies within the limits of the facette mirrors 
    # hexagonal boundaries.
    # If this is the case I have found the mirror, which is hit, and 
    # can calculate:
    # the distance of the intersection point from the center of the facette mirror
    # and the angle relative to the mirror (normal or plane not sure yet)

    # Find intersection of photonline with every mirror-plane
    # calculate intersection point coordinates for every mirror plane
    # mirror plane defined as:
    #   mplane : mplane_n . x = mplane_p . mplane_n
    #   mplane_p is center of mirror
    #   mplane_n is normal vector of mirror
    #
    # photon line defined as:
    #   photon_l0 : x = l0_d * p + l0_p
    #   l0_d is photon direction
    #   l0_p is photon origin
    #   p is free parameter
    #
    #   p0 solves this:
    # --> mplane_n . (l0_d * p + l0_p) = mplane_p . mplane_n
    # --> p * mplane_n . l0_d + mplane_n . l0_p = mplane_p . mplane_n
    # --> p * mplane_n . l0_d  = mplane_n . (mplane_p - l0_p)
    # --> p0 = mplane_n . (mplane_p - l0_p) / (mplane_n . l0_d)
    # 
    # intersection point is called i0 = l0_d *p0 + l0_p
    # but we really need the vector from mi0 = mplane_p to i0 

    logger.debug("   working on keys: %s" % str([pos_key, dir_key]))
    for mindex, mirror in enumerate(mirrors):
        # here I tried to work on a smaller sample of photons, since I 
        # do not have to findMirrorIntersection photons who have already hit.
        # but since photons=photons[...something...] is making a copy
        # or something like that ... no change on photons had any effect :-(
        #print "      current mirror", mindex
        #photons = photons[ np.isnan(photons[:,-1]) ]
        
        #print "working on photons:", run.d['photon_id'].shape[0]
        
        mir = mirror.inner_circle_radius
        mor = mirror.outer_circle_radius
        mplane_n = mirror.dir
        mplane_p = mirror.pos
        
        l0_p = data[pos_key]
        l0_d = data[dir_key]
        
        #~ print "l0_p",l0_p
        #~ print "l0_d",l0_d
        #~ print "l0_p.shape",l0_p.shape
        #~ print "l0_d.shape",l0_d.shape
        #~ 
        #~ print "mirror:",mirror
        #~ print "mplane_n",mplane_n
        #~ print "mplane_p",mplane_p
        #~ print "mplane_n.shape",mplane_n.shape
        #~ print "mplane_p.shape",mplane_p.shape
        #~ 
        
        p0 = -1. * np.dot((l0_p - mplane_p),mplane_n) / np.dot(  l0_d,mplane_n )
        i0 = (l0_d.transpose() * p0).transpose() + l0_p
        
        # the 30. in the following line is simply the speed of light in units of cm/ns
        # so time_needed is in ns of course.
        time_needed = np.sqrt(((i0-l0_p)**2).sum(axis=1)) / 30.
        
        mi0 = i0 - mplane_p
        dist0 = np.sqrt((mi0*mi0).sum(axis=1))
        
        # dist0 lets us decide if:
        #   the point is inside the inner circle of the hexagon? --> Hit
        #   the point is outside the outer circle --> Miss
        #   the point is in beween --> closer look.
        
        # inner circle hits photon_table_indices
        circle_hits_ids = np.where(dist0 <= mirror.inner_circle_radius)[0]
        #print "circle hits: ", len(circle_hits_ids)
        
        #photons[:,-1][dist0 <= mirror.inner_circle_radius] = mindex
        
        # unclear ring hits, where we have to take a closer look
        ring_hits_ids = np.where((dist0 > mir) & (dist0 < mor))[0]
        #print "ring hits: ", len(ring_hits_ids)
        
        mm_mi0 = mi0[ring_hits_ids]
        x_in_plane = np.abs(np.dot(mm_mi0, mirror.x_vector) / np.dot(mirror.x_vector, mirror.x_vector))
        y_in_plane = np.abs(np.dot(mm_mi0, mirror.y_vector) / np.dot(mirror.y_vector, mirror.y_vector))
        
        # these are just the ids in the ids list... crazy huh?
        hex_hits_sub_ids = np.where((x_in_plane < 1) & (y_in_plane < 1. - x_in_plane/2.))[0]
        
        hex_hits_ids = ring_hits_ids[hex_hits_sub_ids]
        #print ring_hits_ids
        #print hex_hits_sub_ids
        #print hex_hits_ids
        #print "hex edge hits: ", len(hex_hits_ids)
        
        #### 
        # ok. now we have all the information, so we store it inside the run
        #
        # we store the mirror ids
        data[new_mirror_id_key][circle_hits_ids] = mindex
        data[new_mirror_id_key][hex_hits_ids] = mindex
        #
        # ... and the coordinates of the intersection points on the mirrors
        data[new_pos_key][circle_hits_ids] = i0[circle_hits_ids]
        data[new_pos_key][hex_hits_ids] = i0[hex_hits_ids]
        #
        # and we store the time needed from the ground to the mirror
        data[new_time_key][circle_hits_ids] = time_needed[circle_hits_ids]
        data[new_time_key][hex_hits_ids] = time_needed[hex_hits_ids]
    
    
    # now find all photons, which have not hit a mirror
    # they still have NaN in ther mirror_id
    # or in other words, find those. which have not isnan
    hit_ids = np.where(~np.isnan(data['mirror_id']))[0]
    
    # now wew delete all those, which have not hit any mirror
    for key in data:
        data[key] = data[key][hit_ids]

    logger.debug("---------- findMirrorIntersection duration: %0.4fs" % (time.time()-start,))
    return run

def reflectOnMirrors( run, dictname='photonData', mirror_alignment_point_key='mirror_alignment_point',
        pos_key='pos_on_mirror', dir_key='ground_dir_turned',
        new_angle_key='angle_on_mirror', new_dir_key='dir_on_mirror' ):
    """  Calculate the new direction of each photon after reflection.
    
    As input one needs:
    
    * the postion of the photon --> pos_key
    * the direction of the photon --> dir_key
    * the direction of the mirror normal --> from mirror_alignment_point_key
    
    From this the new direction is calculated:
    
    * new_dir_key
    
    For convenience the angle to the mirror normal is calculated and stored as well
    
    * new_angle_key

    While *new_dir_key* is the new direction of the photon after its reflection, 
    *new_angle_key* is just the angle of the photon towards the mirror normal, in case somebody wants to know it.
    
    No photons are deleted in this step.
    """
    start = time.time()
    logger.debug("reflectOnMirrors starting up ----------")
    logger.debug("   working on keys: %s" % str([pos_key,dir_key,mirror_alignment_point_key]))
    logger.debug("   Adding new keys to run['%s']: %s" % (dictname,str([new_angle_key, new_dir_key]) ) )
    data = run[dictname]
    mirror_alignment_point = run['settings'][mirror_alignment_point_key]

    # In order to find the new direction of the photon I need:
    #   * current direction &
    #   * normal vector of mirror at current postion of photon
    # 
    # The normal vector I can calculate from the current position
    # Since it is assumed, that every mirror is perfectly aligned.
    # This means for FACTs Davies Cotton reflector, that each mirror is
    # oriented towards a point on the central telescope axis in a distace
    # of twice the mean focal length of all mirrors.
    # Since all mirrors are pherical, their normal vectors at any position 
    # point to that point.
    #
    # Please note, the followings things are not yet taken into account:
    #   * the individual focal length of each mirror.
    #       This would lead to a different point, used for each mirror,
    #       when we construct the normal vector
    #
    #   * the PSF of each mirror.
    #       I do not know, how I should introduce this effect.
    #       I believe, I should calculate the new direction and then add
    #       a little gaussian(?) noise to the result, that's all.
    
    # So let's first calculate the normal vector for 
    # every photons position
    
    i = data[pos_key]
    
    #construct
    n = (i-mirror_alignment_point.pos) * -1.
    # normalize
    n = (n.transpose() / np.sqrt((n**2).sum(axis=1))).transpose()
    
    dir = data[dir_key]
    
    # the angle is not really needed, I just thouht it might be nice to know it.
    n_dot_dir = np.sqrt((dir * n).sum(axis=1))
    angle = np.arccos(n_dot_dir)
    data[new_angle_key] = angle
    # For finding the new direction, I will do as follows:
    #   When you start at the intersection point *i* and follow the 
    #   photon direction one step to *i + dir* you and a new point I call 
    #   "direction point" or *dp*
    #   When you mirror this point at the line throug *i* in *n* direction,
    #   you'll get a new point called "new direction point"  or *ndp*.
    #   The vector from *i* to *ndp* is the new direction vector.
    #
    #   In order to mirror *dp* at that line, I will:
    #       1. construct a plane defined by *dp* and the normal vector *n*
    #       2. construct a line through *i* in *n* direction
    #       3. find crossing point *cp*
    #       4. use vector from *dp* to *cp* to reach *ndp*
    #
    #   dp = i + dir
    #
    #   plane: n . x = n . dp
    #
    #   line: x = n * p + i
    #
    #   find crosspoint:
    #   --> n . (n * p + i) = n . (i + dir)
    #   --> p*n^2 + n.i = n.i + n.dir       // - n.i
    #   --> p * n^2 = n.dir                 // n^2 == 1. since n is normal vector
    #   --> p0  = n.dir
    #   cp = n.dir * n + i
    #
    #   "dp-to-cp" = cp - dp = n.dir * n + i - i + dir
    #              = n.dir * n + dir
    #
    #   ndp = dp + 2 * "dp-to-cp"
    #       = i + 3*dir + 2*n.dir*n
    #
    #   "new direction" = ndp - i = 3*dir + 2*n.dir*n
    new_dir =  dir * -1. + (n.transpose() * n_dot_dir).transpose() * 2.
    
    # the new direction will be saved inside the photons array
    data[new_dir_key] = new_dir
    
    logger.debug("---------- reflectOnMirrors duration: %0.4fs" % (time.time()-start,))
            
    return run



def findFocalPlaneIntersection( run, dictname='photonData', focal_plane_key='focal_plane',
            pos_key='pos_on_mirror', dir_key='dir_on_mirror',
            new_pos_key='pos_in_focal_plane', new_time_key='time_to_focal_plane'):
    """ find new position of photons in focal plane
    
    As input we need information about:
    
    * the focal plane itself --> *focal_plane_key*
    * the current postion of the photon --> *pos_key*
    * and the current direction of the photon --> *dir_key*
    
    As output the new position inside the focal plane is written into *new_pos_key*
    and the time the photon needed is written into *new_time_key*.
    
    The focal plane has a finite size, so all photons, which do not hit the circle 
    are delete in this step.
    """
    
    start = time.time()
    data = run[dictname]
    focal_plane = run['settings'][focal_plane_key]
    logger.debug("findFocalPlaneIntersection starting up ----------")
    logger.debug("   working on keys: %s" % str([pos_key,dir_key,focal_plane_key]))
    logger.debug("   Adding new keys to run['%s']: %s" % (dictname,str([new_pos_key, new_time_key]) ) )
    

    # now the photon has a new direction *new_direction* and is starting
    # from the intersection point *i*
    
    pos = data[pos_key]
    dir = data[dir_key]
    
    fpp = focal_plane.pos
    fpd = focal_plane.dir
    
    p0 = -1. * np.dot((pos - fpp),fpd) / np.dot(  dir,fpd )
    i0 = (dir.transpose() * p0).transpose() + pos
    
    time_needed = np.sqrt(((i0-pos)**2).sum(axis=1)) / 30.
    
    mi0 = i0 - fpp
    dist0 = np.sqrt((mi0*mi0).sum(axis=1))
    
    data[new_pos_key] = i0
    
    # the 30. in the following line is simply the speed of light in units of cm/ns
    # so the time is in ns of course.
    data[new_time_key] = np.sqrt(((i0-pos)**2).sum(axis=1)) / 30.
    
    # Here we delete those, that are not even in a circle with radius=focal_plane.size
    # around the camera center.
    window_hit_ids = np.where ( dist0 < focal_plane.size)[0]
    
    for key in data:
        data[key] = data[key][window_hit_ids]
    logger.debug("---------- findFocalPlaneIntersection duration: %0.4fs" % (time.time()-start,))
    return run

def findPixelChid( run, dictname='photonData', 
        focal_plane_key='focal_plane',
        pos_key='pos_in_focal_plane', 
        new_chid_key='chid', 
        new_Coordinator_key='Coordinator', 
        new_pixel_center_key='pixel_center'):
    """ Sort photons into Pixels accoring to ther current position in the focal plane.
    
    As input we need information about:
    
    * the focal plane itself --> *focal_plane_key*
    * the current postion of the photon --> *pos_key*
    
    As output we have:
    
    * the continous hardware ID (Chid) this photon fell into --> *new_chid_key*
    * the euclidean coordinates of the pixel, with this Chid --> *new_pixel_center_key*
    
    And in addition, we have stored the coor2.Coordinator object in run[new_Coordinator_key],
    since it might be handy to have this object later on.
    A Coordinator knows everything about Pixel IDs and so on.
        
    Giving it the focal plane is needed, since the *pos_in_focal_plane* is stored in real 3D coordinates
    and for sorting the photons into the pixels I need their positions in 2D coordinates inside the plane.
    In this special case I simply could drop the z-coordinate, but in more general cases I would need the focal_plane.
    Maybe this shows me, that after the focal_plane was needed for the first time, it should be introduced into
    the data flow, so that every upcoming processor can make use of it.

    After this call we have deleted all those photons, which did not fall into a Pixel.
    """
    
    from coor2 import Coordinator
    start = time.time()
    data = run[dictname]
    focal_plane = run['settings'][focal_plane_key]

    logger.debug("findPixelChid starting up ----------")
    logger.debug("   working on keys: %s" % str([pos_key,focal_plane_key]))
    logger.debug("   Adding new keys to run['%s']: %s" % (dictname,str([new_chid_key, new_Coordinator_key, new_pixel_center_key]) ) )

    # calculate the position in the focal plane relative to the 
    # center of the focal plane center. 
    pos = data[pos_key][:] - focal_plane.pos
    # After that step, if the telescope is not turned all 
    # postions should have zero as last coordinate
    # so lets throw an Exception if this is not the case.
    np.testing.assert_almost_equal(pos[:,-1], np.zeros(len(pos[:,-1])))
    
    # I only want x-y-coordinates, not x-y-z
    pos = pos[:,:-1]
    
    c = Coordinator()
    chids = c.get_chid_from_pos(pos)

    data[new_chid_key] = chids[:]
    
    # lets also store the coordinator inside the run somewhere
    run[new_Coordinator_key] = c
    
    # and finally lets delete all those, that did not fall into a hexagon.
    fell_in_hexagons = np.where( ~np.isnan(chids) )[0]
    for key in data:
        data[key] = data[key][fell_in_hexagons]
    
    indices_for_coordinator = c.chid_to_array_index[data[new_chid_key].astype(int)]
    data[new_pixel_center_key] = c.pos[indices_for_coordinator]
    
    logger.debug("---------- findPixelChid duration: %0.4fs" % (time.time()-start,))
    return run


def _read_winston_cone_file ( filename ):
    """ the file looks like this:
    
        #   softID:  0 - 1439  (spiral from center)
        #   hardID:  CBPX    Crate 0..3  Board 0..9  Patch 0..3  piXel 0..8
        #   G-APD:   G_APD chip identifier
        #   V_op:    Bias Voltage for G-APD (including Bias offset)
        #   HV_B:    Bias Voltage Board   0..9
        #   HV_Ch:   Bias Voltage Channel 0..32
        #   pos_X/pos_Y: position in pixel-size units
        #
        #           CBPX
        # softID  hardID  geom_i  geom_j   G-APD   V_op  HV_B HV_C:  pos_X     pos_Y     radius
        0    1036       0       0     957   70.59   4   7      0.00      0.50      0.25
        1    1033       0       1     997   70.62   4   6      0.00     -0.50      0.25
        2    1034      -1       0    1087   70.59   4   7     -0.87      0.00      0.75
    """
    
    pixel = []
    f = open( filename )
    
    # okay, so the co called hardID is CBPX .. but I need the so called CHid
    # which I get by:
    # C*360 + B*36 + P*9 + X
    
    chid_to_hex_coors = {}
    
    # lets form the pos in a affine coordinate system
    pixel_size = 0.9 # cm
    
    center = np.array([0, pixel_size/2., 0.])
    e1 = np.array([0, pixel_size, 0.])
    e2 = np.array([np.sqrt(3./4.)*pixel_size, pixel_size/2., 0.])
    
    euclidic_coors = np.zeros( (1440,3), dtype=np.float32)
    hex_coors = np.zeros( (1440,2), dtype=np.int32)
    chid = np.zeros( 1440, dtype=np.int32)
    
    for line in f:
        line = line.split('#')[0]
        line = line.split()
        if not line:
            continue
        
        soft_id = int(line[0]) # can be used as line counter :-)
        
        CBPX = int(line[1])
        X = CBPX % 10
        P = (CBPX - X) % 100 / 10
        B = (CBPX - X - 10*P) % 1000 / 100
        C = (CBPX - X - 10*P - 100*B) / 1000
        chid_i = C*360 + B*36 + P*9 + X
        
        chid[soft_id] = chid_i
        euclidic_coors[soft_id,0] = float(line[8]) * pixel_size
        euclidic_coors[soft_id,1] = float(line[9]) * pixel_size
        euclidic_coors[soft_id,2] = 0.
        
        ec = euclidic_coors[soft_id]
        #hex_coors[soft_id,0] = np.round(np.dot(ec-center, e1)).astype(np.int32)
        #hex_coors[soft_id,1] = np.round(np.dot(ec-center, e2)).astype(np.int32)
        hex_coors[soft_id,0] = np.round(2*np.dot(ec-center, e1) / np.linalg.norm(e1)**2)
        hex_coors[soft_id,1] = np.round(2*np.dot(ec-center, e2) / np.linalg.norm(e2)**2)
    
    #######################
    #
    # create look up dict
    hex_coors_to_chid = {}
    for index, hc in enumerate(hex_coors):
        key = tuple(hc.tolist())
        hex_coors_to_chid[key] = chid[index]
    
    all = {}
    all['euclidic_coors']  =euclidic_coors
    all['hex_coors']  =hex_coors
    all['chid']  =chid
    all['hex_coors_to_chid']=hex_coors_to_chid
    
    
    return all



def _reflect_photon( photon, mirrors):
    """ finds out: 
            which mirror is hit by photon 
            and where 
            and in which angle relative to mirror 
    """
    
    # the line defined by the photon is used to find the intersection point 
    # with the plane of each facette mirror. Then I check,
    # if the intersection point lies within the limits of the facette mirrors 
    # hexagonal boundaries.
    # If this is the case I have found the mirror, which is hit, and 
    # can calculate:
    # the distance of the intersection point from the center of the facette mirror
    # and the angle relative to the mirror (normal or plane not sure yet)
    
    for mirror in mirrors:
        #facette mirror plane, defined as n . x = d1 . n
        n = mirror.dir
        d1 = mirror.pos
        
        # line of photon defined as r = lambda * v + d2
        v = photon.dir
        d2 = photon.pos
        
        # the intersection coordinates are found by solving
        # n . (lambda * v + d2) - d1 . n == 0, for lambda=lambda_0
        # and then the intersection is: i = lambda_0 * v + d2
        #
        # putting int in another form:
        # solve:
        # lambda * n.v + n.d2 - n.d1 == 0
        # or
        # lambda_0 = n.(d1-d2) / n.v
        
        # FIXME: if one of the two dot-products is very small,
        # we shuold take special care maybe
        # if n.(d1-d2) is very small, this means that the starting point of 
        #   the photon is already nearly in the plane, so lambda_0 is expected to
        #   be very small ... erm .. maybe this is actually not a special case
        #   but very good.
        # of n.v is very small, this means the patch of the photon is nearly
        #   parallel to the plane, so the error ob lambda_0 might be very large.
        #   in addition, this might just tell us, that the mirror is hit under 
        #   strange circumstances ... so its not a good candidate and we can already go on.
        lambda_0 = (np.dot(n,(d1-d2)) / np.dot(n,v))
        
        #intersection between line and plane
        i = lambda_0 * v + d2
        
        # I want the distance beween i and d1 so I can already from the distance find 
        # out if this is our candidate.
        i_inside_plane = i-d1
        distance = helpers.length(i_inside_plane)

        # now lets check if the mirror intersection point really lies inside 
        # the hexagon or not.
        # first we do two simple test on the distance.
        # if i inside the hexagons inner circle, we are positive already.
        # no further tests needed
        if distance <= mirror.inner_circle_radius:
            break
        # if on the other hand i is way too far out, we are clearly negative
        elif distance >= mirror.outer_circle_radius:
            mirror = None
        # in all other cases we have to look closer.
        # so we will calculate the points "x" and "y" coordinates inside the 
        # hexagon plane. The x-axis should look eastwards, when the telescope 
        # look up to zenith. and the hexagon y-axis should look northwards.
        # since the mirrors are turned, we can not simply take the normal
        # x and y-axis vectors, but each mirror has his own normal vectors.
        # also these vectors are not of length one, but:
        # the hexagon x-axis vector goes from center to the middle of the right edge
        # and the y-axis vector goes from the center to the upper corner.
        else:
            x_in_plane = np.dot(i_inside_plane, mirror.x_vector) / np.dot(mirror.x_vector, mirror.x_vector)
            y_in_plane = np.dot(i_inside_plane, mirror.y_vector) / np.dot(mirror.y_vector, mirror.y_vector)
            if np.abs(x_in_plane) > 1:
                mirror = None
            elif np.abs(y_in_plane) > (1. - np.abs(x_in_plane)/2.):
                mirror = None
            else:
                # we are inside the hexagon, so we have found the mirror we were searching for
                break
    if not mirror is None:
        photon.mirror_index = mirror.index
        photon.mirror_intersection = i
        photon.mirror_center_distance = distance
        #print "mirror found:", mirror.index , 
        #print "distance", distance
        # now I have to find out, if the photon is not only in the 
        # right distance but actually has hit the mirror.
        # this I do like this
        # i-d1    is a vector in the mirror plane pointing from d1 to the intersection point i.
        # if I know turn the entire mirror plane so it lies withing the x-y-plane
        # by applying a simple turning-matrix, then each vector inside the plane with turn into
        # a nice x,y vector. 
        # now I assume, that the hexagon is "pointing" lets say to into y direction
        # so I can e.g. say:
        # x has to be between -30.3 and +30.3 and y has to be
        # between 35 - m * |x| and -35 + m * |x| ... pretty simple.
        # maybe one can leave the turning aside, but I like that I can imagine it very nicely
        # 
        #
        # I don't do this yet .. since I don't know by heart how a turning matrix looks :-)
        # so I just simulate round mirrors
        ######################################################################
        
        
        # next step, since I know the intersection point, is the new direction.
        # So I need the normal of the mirror in the intersection point.
        # since the normal of every mirror is alway pointing to the camera center
        # this is not difficult.
        
        normal_at_intersection = (mirror_alignment_point.pos - i) / helpers.length(mirror_alignment_point.pos - i)
        #print "normal_at_intersection",normal_at_intersection
        
        angle = np.arccos(np.dot( v, normal_at_intersection) / (helpers.length(v) * helpers.length(normal_at_intersection)))
        photon.angle_to_mirror_normal = angle
        #print "angle:", angle/np.pi*180., "deg"
        
        
        # okay, now I have the intersection *i*, 
        # the old direction of the photon *v*
        # and the normalvector at the intersection.
        ######################################################################
        ######################################################################
        # I do this now differently.
        # I will mirror the "point" at the tip of *v* at the line created by
        # the normalvector at the intersection and the intersection.
        # this will gibe me a mirrored_point *mp* and the vector from *i* to *mp*
        # is the *new_direction* it should even be normalized.
        
        # 1. step: create plane through the "tip" of *v* and the normal_at_intersection.
        # 2. step: find crossingpoint on line through *i* and the normal_at_intersection,
        # 3. step: vector from "tip" of *v* to crossingpoint times 2 points to 
        #           the "tip" of *mirrored_v*
        
        # plane: n_plane_3 . r = p_plane_3 . n_plane_3
        # p_plane_3 = i+v
        # n_plane_3 = normal_at_intersection
        
        # line: r = lambda_3 * v_line_3 + p_line_3
        # p_line_3 = i
        # v_line_3 = normal_at_intersection
        
        # create crossing: n_plane_3 . (lambda_3 * v_line_3 + p_line_3) = p_plane_3 . n_plane_3
        #   <=> lambda_3 = (p_plane_3 - p_line_3 ).n_plane_3  / n_plane_3 . v_line_3
        #   <=> lambda_3 = (i+v - i).normal_at_intersection  / normal_at_intersection . normal_at_intersection
        #   <=> lambda_3 = v.normal_at_intersection
        
        lambda_3 = np.dot(v, normal_at_intersection)
        #print "lambda_3", lambda_3
        crossing_point_3 = lambda_3 * normal_at_intersection + i
        #print "crossing_point_3", crossing_point_3
        
        from_tip_of_v_to_crossing_point_3 = crossing_point_3 - (i+v)
        
        tip_of_mirrored_v = i+v+ 2*from_tip_of_v_to_crossing_point_3
        
        new_direction = tip_of_mirrored_v - i
        
        #print "new_direction",new_direction
        #print "old direction", v
        photon.new_direction = new_direction
        ######################################################################
        ######################################################################
        """
        # both directions form a plane, and when I turn the old *v* by
        # twice the angle between *v* and *normal_at_intersection* 
        # inside this plane then I get the new direction of the photon.
        
        # so lets first get the normal of the reflection plane
        normal_of_reflection_plane =np.cross( v, normal_at_intersection)
        
        print length(normal_of_reflection_plane), "should be one"
        print length(v), "should be one"
        print length(normal_at_intersection), "should be one"
        print np.dot(v, normal_at_intersection), "should *NOT* be zero"
        print np.dot(v, normal_of_reflection_plane), "should be zero"
        print np.dot(normal_at_intersection, normal_of_reflection_plane), "should be zero"
        
        angle = np.arccos(np.dot( v, normal_at_intersection) / (length(v) * length(normal_at_intersection)))
        photon.angle_to_mirror_normal = angle
        print "angle:", angle/np.pi*180., "deg"
        
        # the rotation matrix for the rotation of *v* around normal_of_reflection_plane is
        R = make_rotation_matrix( normal_of_reflection_plane, 2*angle )
        
        print "R"
        pprint(R)
        
        new_direction = np.dot( R, v)
        photon.new_direction = new_direction
        
        print "old direction", v, length(v)
        print "new direction", new_direction, length(new_direction)
        print "mirror center", mirror.pos
        print "interception point", i
        print "center of focal plane", focal_plane.pos
        """

        # now the photon has a new direction *new_direction* and is starting
        # from the intersection point *i*
        # now I want to find out where there focal plane is hit.
        # So I have to repeat the stuff from up there
        
        #print "np.dot(focal_plane.dir,new_direction))", np.dot(focal_plane.dir,new_direction)
        
        lambda_1 = (np.dot(focal_plane.dir ,(focal_plane.pos - i)) / np.dot(focal_plane.dir,new_direction))
        
        #print "lambda_1", lambda_1
        focal_plane_spot = lambda_1 * new_direction + i
        #print "focal_plane_spot",focal_plane_spot
        photon.hit = True
        focal_plane_pos = focal_plane_spot - focal_plane.pos
        photon.focal_plane_pos =focal_plane_pos
        #photon.hit = True
        if helpers.length(focal_plane_pos) <= focal_plane.size:
            photon.hit = True
        else:
            photon.hit = False
            return photon
        # now as a final step we have to find the coordinates of the vector
        # from the center of the focal plane to the spot where the photon 
        # actually hit the focal plane, as if the plane was not turned.
        # so if we turn the plane back into the x-y-plane
        # our *focal_plane_pos* vector has only two coordinates x,y, which are non-zero.
        # so lets do that.
        # in order to do so, we need the angles, by which the telescope
        # was turned .. we hae made them global variables !!ugly i know!!
        
        R = make_rotation_matrix( np.array([0,1,0]), -1.*telescope_theta/ 180. *np.pi )
        turned_focal_plane_pos = np.dot( R, focal_plane_pos)
        R = make_rotation_matrix( np.array([-1,0,0]), -1.*telescope_phi/ 180. *np.pi )
        turned_focal_plane_pos = np.dot( R, turned_focal_plane_pos)
        photon.turned_focal_plane_pos = turned_focal_plane_pos
        #if np.abs(turned_focal_plane_pos[2] ) > 1e-12:
        #    print turned_focal_plane_pos[2]
        #    raise Exception("the z-coordinate should be zero but is larger than 1e-12")
        
        
        
        #print "distance from focal plane center=",  length(focal_plane_spot-focal_plane.pos)
    else:
        photon.hit = False
    return photon
    

def _sort_into_pixel( run, pixel ):
    """
    """
    
    # Inside the camera plane one can define a Coordinate System like this
    pixel_size = 0.95 # cm
    center = np.array([0., 0.5])*pixel_size
    ey = np.array([0., -1.])*pixel_size
    ex = np.array([np.sqrt(3)/2. , 0.5])*pixel_size
    # the center of each Pixel has integer coordinates in this coordinate system.
    # If I now transfor the coordinates of all photons, I will get in most cases
    # non integer values.
    # Now it would be wrong to assume, that I may just round the coordinates
    # to the next integer, in order to find the pixel, where this photon 
    # fell into.
    # But I can already find 4 pixels, where the one I search is in.
    #
    # So lets first calculate the photons coordinates in this hexagonal coordinate system
    
    # calculate distance of each pixel to each photon point
    pos = run.d['pos_in_focal_plane'][:] - focal_plane.pos
    # I only want x-y-coordinates, not x-y-z
    pos = pos[:,:-1]
    
    # Now we do the transformation:
    # FIXME: there must be a nicer more general way to write this down, but
    # I am to stupid. 
    #
    hexAffineCoor = pos.copy()
    pos_shifted = pos - center
    hexAffineCoor[:,0] = pos_shifted[:,0] / ex[0]
    pos_shifted -= np.outer(hexAffineCoor[:,0],ex)
    hexAffineCoor[:,1] = pos_shifted[:,1] / ey[1]
    #
    # Ok now we have the hexagonal coordinates of each photon
    # The coordinates of the 4 candiate pixels are now generated
    # by rounding the hexagonal coordinates up and down... 
    # or calculating the floor and ceil ... however you want to call it
    cand_x_min_hex = np.floor(hexAffineCoor[:,0]).astype(int)
    cand_x_max_hex = np.ceil(hexAffineCoor[:,0]).astype(int)
    cand_y_min_hex = np.floor(hexAffineCoor[:,1]).astype(int)
    cand_y_max_hex = np.ceil(hexAffineCoor[:,1]).astype(int)
    #
    # okay so the 4 candidates for each pixel are:
    # cand1_hex = (cand_x_min , cand_y_min)
    # cand2_hex = (cand_x_min , cand_y_max)
    # cand3_hex = (cand_x_max , cand_y_min)
    # cand4_hex = (cand_x_max , cand_y_max)
    #
    # Now I calculate for each of them the normal euclidian coordinates again
    cand1_euc = np.outer(cand_x_min,ex) + np.outer(cand_y_min,ey) + center
    cand2_euc = np.outer(cand_x_min,ex) + np.outer(cand_y_max,ey) + center
    cand3_euc = np.outer(cand_x_max,ex) + np.outer(cand_y_min,ey) + center
    cand4_euc = np.outer(cand_x_max,ex) + np.outer(cand_y_max,ey) + center
    #
    #
    # Okay and finally I will calculate the distance of the pixel for each candidate:
    # Such that I get a matrix like: num_pixel x num_candidates
    # And when I seach for the 
    
    
    
    
    
    
    
    
    """
    run.d['nearest_pixel'] = np.zeros( pos.shape[0], dtype=np.int32)
    run.d['dist_to_pixel'] = np.zeros( pos.shape[0], dtype=np.int32)
    
    # a loop over all photons
    for index,p in enumerate(pos):
        distance_vectors = pixel - p
        distances = np.sqrt((distance_vectors**2).sum(axis=1))
        nearest_pixel = distances.argmin()
        dist_to_pixel = distances.min()
        run.d['nearest_pixel'][index] = nearest_pixel
        run.d['dist_to_pixel'][index] = dist_to_pixel
        
    
    
    # lets form the pos in a affine coordinate system
    #pixel_size = 0.9 # cm
    #center = np.array([0, pixel_size/2., 0.])
    #e1 = np.array([0, pixel_size, 0.])
    #e2 = np.array([np.sqrt(3./4.)*pixel_size, pixel_size/2., 0.])
    
    #pos = pos - center
    #coor1 = np.round(np.dot(pos, e1)).astype(int)
    #coor2 = np.round(np.dot(pos, e2)).astype(int)
    #return coor1,coor2
    """
    return run



def reflector( run ):
    """ Calculate the pixel ID and the time for each photon.
    
    
    """
    
    run = pushMirrorInfo( run )

    run = findMirrorIntersection( run, dictname='photonData', 
            pos_key='ground_pos_turned', 
            dir_key='ground_dir_turned',
            mirror_key='mirrors',
            new_pos_key='pos_on_mirror', 
            new_mirror_id_key='mirror_id', 
            new_time_key='time_to_mirror')

    run = reflectOnMirrors( run, dictname='photonData', 
            mirror_alignment_point_key='mirror_alignment_point',
            pos_key='pos_on_mirror', 
            dir_key='ground_dir_turned',
            new_angle_key='angle_on_mirror', 
            new_dir_key='dir_on_mirror' )
        
    run = findFocalPlaneIntersection( run, dictname='photonData', 
            focal_plane_key='focal_plane',
            pos_key='pos_on_mirror', 
            dir_key='dir_on_mirror',
            new_pos_key='pos_in_focal_plane', 
            new_time_key='time_to_focal_plane')
            
    run = findPixelChid( run, dictname='photonData', 
            focal_plane_key='focal_plane',
            pos_key='pos_in_focal_plane', 
            new_chid_key='chid', 
            new_Coordinator_key='Coordinator', 
            new_pixel_center_key='pixel_center')

    #pixelstuff = read_winston_cone_file ( "../tel_definitions/FACTmap111030.txt" )
    #sort_into_pixel( run, pixelstuff['euclidic_coors'] )
    #
    #coor = Coordinator("map.txt")
    
    #c1,c2 = sort_into_pixel( run, pixel )
    
    #run = main( run )
    return run

def make_plots( run,dictname='photonData' ):
    data = run[dictname]
    helpers.plot_in_TH2F( data['ground_pos'][:,:-1], 'ground_pos' )
    helpers.plot_in_TH2F( data['pos_on_mirror'][:,:-1], 'pos_on_mirror' )
    helpers.plot_in_TGraph2D( data['pos_on_mirror'] , 'pos_on_mirror as plot')
    helpers.plot_in_TGraph2D( data['ground_dir'] ,'ground_dir as plot')
    helpers.plot_in_TGraph2D( data['ground_dir_turned'] ,'ground_dir_turned')
    helpers.plot_in_TGraph2D( data['dir_on_mirror'] ,'dir_on_mirro')
    helpers.plot_in_TH1F( data['angle_on_mirror'] ,'angle_on_mirror')
    helpers.plot_in_TH2F( data['dir_on_mirror'][:,:-1] ,'dir_on_mirror')
    helpers.plot_in_TH2F( data['pos_in_focal_plane'][:,:-1] ,'pos_in_focal_plane', bins=(100,100) )
    #helpers.plot_in_TH2F( data['pixel_center'] ,'pixel_center',bins=(50,50))
    
    # This is a module from PyFACT
    import plotters
    camplotter = plotters.CamPlotter("Pixel Center", dpi=80, s=80)
    
    photons_in_pixel = np.zeros(1440)
    for i in range(1440):
        photons_in_pixel[i] = (data['chid'] == i).sum()
        
    camplotter(photons_in_pixel)
    plt.xlim(-22,22)
    plt.ylim(-22,22)
    
    

if __name__ == '__main__':
    import rlcompleter
    import readline
    readline.parse_and_bind('tab: complete')
    from optparse import OptionParser

    parser = OptionParser()
    parser.add_option("-f", "--file", dest="filename",
        help="write report to FILE", metavar="FILE")
    parser.add_option("-p", "--plot",
        action="store_true", dest="create_plots", default=False,
        help="create some informative debugging plots")
    options, args = parser.parse_args()
    if options.filename is None:
        print 'Please submit a filename. "', sys.argv[0],'-h" for help'
        sys.exit(1)

    import readcorsika
    import absorb_photons as ap    
    run = readcorsika.readCorsikaFile( options.filename )
    run = ap.absorb_photons(run)
    
    run = reflector(run)
    
    if options.create_plots:
        make_plots(run)
